MODULE SRFSN_WEBALSTL_MOD
CONTAINS
SUBROUTINE SRFSN_WEBALSTL(KIDIA,KFDIA,KLON,KLEVSN,&
 & PTMST  , PFRTI     , PCIL     , LDLAND, LLNOSNOW, &
 & PSSNM1M5, PWSNM1M5 , PRSNM1M5 , PTSNM1M5,&
 & PSLRFL5 , PSSRFLTI5, PAHFSTI5 , PEVAPTI5,&
 & PSSFC5  , PSSFL5   , PEVAPSNW5, &
 & PTSFC5  , PTSFL5   ,           &
 & PTSURF5 , PTICE5, PSNOTRS5 ,           &
 & PAPRS5  , PWSURF5  , &
 & PSSDP3,&
 & YDSOIL , YDCST   ,             &
 & PTSN5, PGSN5,                  &
!--perturbations
 & PSSNM1M,PWSNM1M,PTSNM1M,                 &
 & PSLRFL,PSSRFLTI,PAHFSTI ,PEVAPTI,        &
 & PSSFC, PSSFL, PEVAPSNW, PTSFC  ,PTSFL   ,&
 & PTSURF, PTICE, PSNOTRS,                         &
 & PAPRS, PWSURF,                           &
 & PTSN, PGSN)


USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK

USE YOS_SOIL , ONLY : TSOIL 
USE YOS_CST  , ONLY : TCST
USE EC_LUN   , ONLY : NULERR

USE ABORT_SURF_MOD
USE YOMSURF_SSDP_MOD

#ifdef DOC
! (C) Copyright 2021- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!**** *SRFSN_WEBALSTL* - CONTAINS SNOW PARAMETRIZATION
!                        (Tangent linear)
!
!     PURPOSE.
!     --------
!          COMPUTES CHANGES IN SNOW TEMPERATURE

!**   INTERFACE.
!     ----------
!          *SRFSN_WEBALSTL* IS CALLED FROM *SURFTSTPSTL*.

!     PARAMETER   DESCRIPTION                                          UNITS
!     ---------   -----------                                          -----

!     INPUT PARAMETERS (INTEGER):
!     *KIDIA*      START POINT
!     *KFDIA*      END POINT
!     *KLON*       NUMBER OF GRID POINTS PER PACKET

!     INPUT PARAMETERS (REAL):
!     *PTMST*      TIME STEP                                            S
!     *PFRTI*      TILE FRACTIONS                                       -
!     *PCIL*       FRACTION OF LAND GLACIER COVERED
 
!     INPUT PARAMETERS (LOGICAL):
!     *LDLAND*      TRUE WHEN LAND POINT
!     *LLNOSNOW*    TRUE WHEN SNOW IS NOT PRESENT

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PSSNM1M5    PSSNM1M       SNOW MASS (per unit area)                  kg/m2/s
!  PTSNM1M5    PTSNM1M       SNOW TEMPERATURE                           K
!  PRSNM1M5     ---           SNOW DENSITY                               kg/m3 
!  PDSNM1M5    PDSNM1M       SNOW TOTAL DEPTH (in m)                    m
!  PTSAM1M5    PTSAM1M       SOIL TEMPERATURE                           K
!  PHLICEM1M5  ---           LAKE ICE THICKNESS                         m
!  PSLRFL5     PSLRFL        NET LONGWAVE  RADIATION AT THE SURFACE     W/m2
!  PSSRFLTI5   PSSRFLTI      NET SHORTWAVE RADIATION AT THE SURFACE,    W/m2
!                            FOR EACH TILE  
!  PAHFSTI5    PAHFSTI       TILE SENSIBLE HEAT FLUX                    W/m2
!  PEVAPTI5    PEVAPTI       TILE EVAPORATION                           kg/m2/s
!  PSSFC5      PSSFC         CONVECTIVE  SNOW FLUX AT THE SURFACE       kg/m2/s
!  PSSFL5      PSSFL         LARGE SCALE SNOW FLUX AT THE SURFACE       kg/m2/s
!  PEVAPSNW5   PEVAPSNW      EVAPORATION FROM SNOW UNDER FOREST         kg/m2/s
!  PTSFC5      PTSFC         Convective Throughfall at the surface      kg/m2/s
!  PTSFL5      PTSFL         Large Scale Throughfall at the surface     kg/m2/s

!     PARAMETERS AT T+1 :
!  Trajectory  Perturbation  Description                                Unit
!  PTSN5       PTSN          SNOW TEMPERATURE                           K

!    FLUXES FROM SNOW SCHEME:
!  Trajectory  Perturbation  Description                                Unit
!  PGSN5       PGSN          GROUND HEAT FLUX FROM SNOW DECK TO SOIL    W/m2 (#)


! (#) THOSE TWO QUANTITIES REPRESENT THE WHOLE GRID-BOX. IN RELATION
!       TO THE DOCUMENTATION, THEY ARE PGSN=Fr_s*G_s

!     METHOD.
!     -------
!     Adaptation of srfsn_webals (tangent linear)
!     - Initial Liquid water set to zero so far
!     - Interception of rainfall

!     EXTERNALS.
!     ----------

!     REFERENCE.
!     ----------
!          SEE SOIL PROCESSES' PART OF THE MODEL'S DOCUMENTATION FOR
!     DETAILS ABOUT THE MATHEMATICS OF THIS ROUTINE.

!     Original:
!     ---------
!     G. Arduini             E.C.M.W.F.     01-07-2020  

!     Modifications
!     -------------
!     I. Ayan-Miguez (BSC) Sep 2023 Added PSSDP3 object for spatially distributed parameters 
!     ------------------------------------------------------------------
#endif

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON
INTEGER(KIND=JPIM), INTENT(IN)   :: KLEVSN
REAL(KIND=JPRB)   , INTENT(IN)   :: PTMST
REAL(KIND=JPRB),    INTENT(IN)   :: PFRTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCIL(:)

LOGICAL           , INTENT(IN)   :: LDLAND(:)
LOGICAL           , INTENT(IN)   :: LLNOSNOW(:)

REAL(KIND=JPRB)   , INTENT(IN)   :: PSSNM1M5(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSNM1M5(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PRSNM1M5(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSNM1M5(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSURF5(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTICE5(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PSNOTRS5(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PAPRS5(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSURF5(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PSLRFL5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSRFLTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PAHFSTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFC5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFL5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPSNW5(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFC5(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFL5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFC5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFL5(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PSSDP3(:,:,:)

TYPE(TSOIL)       , INTENT(IN)   :: YDSOIL
TYPE(TCST)        , INTENT(IN)   :: YDCST

REAL(KIND=JPRB)   , INTENT(OUT)  :: PTSN5(:,:)
REAL(KIND=JPRB)   , INTENT(OUT)  :: PGSN5(:)


! Perturbation variables
REAL(KIND=JPRB)   , INTENT(IN)   :: PSSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSURF(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTICE(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PSNOTRS(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PAPRS(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSURF(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PSLRFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSRFLTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PAHFSTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFC(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPSNW(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFC(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFL(:)
REAL(KIND=JPRB),    INTENT(IN):: PTSFC(:)
REAL(KIND=JPRB),    INTENT(IN):: PTSFL(:)

REAL(KIND=JPRB)   , INTENT(OUT)  :: PTSN(:,:)
REAL(KIND=JPRB)   , INTENT(OUT)  :: PGSN(:)

! Local variables (in srfsn_driver in full non linear)
REAL(KIND=JPRB) :: ZSNOWF5(KLON)
REAL(KIND=JPRB) :: ZRAINF5(KLON)
REAL(KIND=JPRB) :: ZEVAPSN5(KLON)
REAL(KIND=JPRB) :: ZHFLUX5(KLON)
REAL(KIND=JPRB) :: ZFRSN5(KLON)
REAL(KIND=JPRB) :: ZMELTSN5(KLEVSN)
REAL(KIND=JPRB) :: ZFREZSN5(KLEVSN)

REAL(KIND=JPRB) :: ZTBOTTOM5(KLON)

REAL(KIND=JPRB) :: ZSNOWF(KLON)
REAL(KIND=JPRB) :: ZRAINF(KLON)
REAL(KIND=JPRB) :: ZEVAPSN(KLON)
REAL(KIND=JPRB) :: ZHFLUX(KLON)
REAL(KIND=JPRB) :: ZFRSN(KLON)
REAL(KIND=JPRB) :: ZMELTSN(KLEVSN)
REAL(KIND=JPRB) :: ZFREZSN(KLEVSN)

REAL(KIND=JPRB) :: ZTBOTTOM(KLON)

! Local variables      traj                pert
REAL(KIND=JPRB) :: ZDSN5(KLEVSN)      , ZDSN(KLEVSN)    ! actual snow depth
REAL(KIND=JPRB) :: ZSNHC5(KLEVSN)     , ZSNHC(KLEVSN)   ! snow heat capacity 
REAL(KIND=JPRB) :: ZICE5(KLEVSN)      , ZICE(KLEVSN)    ! snow ice content (PSSN-PWSN)
REAL(KIND=JPRB) :: ZSNCOND5(KLEVSN)   , ZSNCOND(KLEVSN) ! Snow thermal conductivity 
REAL(KIND=JPRB) :: ZSNCONDH5(KLEVSN+1), ZSNCONDH(KLEVSN+1) ! THERMAL CONDUCTIVITY IN HALF LEVEL TERM 
REAL(KIND=JPRB) :: ZTSTAR5(KLEVSN)    , ZTSTAR(KLEVSN)  ! New snow temperature 
REAL(KIND=JPRB) :: ZISTAR5(KLEVSN)    , ZISTAR(KLEVSN)  ! New ice content 
REAL(KIND=JPRB) :: ZWSTAR5(KLEVSN)    , ZWSTAR(KLEVSN)  ! New liquid water content 
REAL(KIND=JPRB) :: ZLIQF5(0:KLEVSN)   , ZLIQF(0:KLEVSN) ! LIQUID WATER FLUX 

REAL(KIND=JPRB) :: ZTA5(KLEVSN),ZTB5(KLEVSN),ZTC5(KLEVSN),ZTR5(KLEVSN)  ! TERMS TO TRI-DIAG
REAL(KIND=JPRB) :: ZTA(KLEVSN), ZTB(KLEVSN),ZTC(KLEVSN),ZTR(KLEVSN)  ! TERMS TO TRI-DIAG
REAL(KIND=JPRB) :: ZGSN5,ZWCAP5,  ZSNVCOND5, ZTMP15
REAL(KIND=JPRB) :: ZGSN, ZWCAP,   ZSNVCOND,  ZTMP1
REAL(KIND=JPRB) :: ZGSNRES5(0:KLEVSN),ZGSNRES(0:KLEVSN) 
REAL(KIND=JPRB) :: ZQ5(1:KLEVSN), ZQ(1:KLEVSN), ZTMP05(1:KLEVSN), ZTMP0(1:KLEVSN)

REAL(KIND=JPRB) :: ZFF5, ZWU5, ZLWT5,  ZINVWSAT5, ZLIC5, ZLAMBDASAT5, ZCOND5, ZKERSTEN5, ZSOILCOND5, ZSURFCOND5
REAL(KIND=JPRB) :: ZFF,  ZWU,  ZLWT,   ZINVWSAT,  ZLIC,  ZLAMBDASAT,  ZCOND,  ZKERSTEN,  ZSOILCOND, ZSURFCOND

REAL(KIND=JPRB) :: ZTMST,ZIHCAP,ZSOILDEPTH1,ZEPSILON
REAL(KIND=JPRB) :: ZWHCAP

REAL(KIND=JPRB) :: ZTMPHC, ZTMPCOND, ZSNCONDH1, ZSNCONDH2, ZCHANGE5,ZCONSLWC, ZLN10

INTEGER(KIND=JPIM) :: JL,JK,KLACT

REAL(KIND=JPRB) :: ZZLIQF5

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

#include "fcsurf.h"

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALSTL_MOD:SRFSN_WEBALSTL',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------
ASSOCIATE(RHOCI=>YDSOIL%RHOCI,RHOICE=>YDSOIL%RHOICE,&
 & RTT=>YDCST%RTT,RLMLT=>YDCST%RLMLT,&
 & RALAMSN=>YDSOIL%RALAMSN, RLAMICE=>YDSOIL%RLAMICE, &
 & RFRTINY=>YDSOIL%RFRTINY, RFRSMALL=>YDSOIL%RFRSMALL, &
 & RDSNMAX=>YDSOIL%RDSNMAX, SNHCONDAV=>YDSOIL%SNHCONDAV, &
 & SNHCONDBV=>YDSOIL%SNHCONDBV, SNHCONDCV=>YDSOIL%SNHCONDCV, &
 & SNHCONDPOV=>YDSOIL%SNHCONDPOV, &
 & RKERST1=>YDSOIL%RKERST1, RKERST2=>YDSOIL%RKERST2, RKERST3=>YDSOIL%RKERST3, &
 & RLAMBDADRY=>YDSOIL%RLAMBDADRY, RLAMBDADRYM3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMBDADRYM3D), &
 & RLAMBDAICE=>YDSOIL%RLAMBDAICE, RLAMBDAWAT=>YDSOIL%RLAMBDAWAT, &
 & RLAMSAT1=>YDSOIL%RLAMSAT1, RLAMSAT1M3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMSAT1M3D), RSIMP=>YDSOIL%RSIMP, &
 & RWSAT=>YDSOIL%RWSAT, RWSATM3D=>PSSDP3(:,:,SSDP3D_ID%NRWSATM3D), &
 & RCONDSICE=>YDSOIL%RCONDSICE, &
 & RLWCSWEA=>YDSOIL%RLWCSWEA , RLWCSWEB=>YDSOIL%RLWCSWEB, RLWCSWEC=>YDSOIL%RLWCSWEC)
! RLMLT (latent heat of fusion J Kg -1)
ZTMST      = 1.0_JPRB/PTMST 
ZIHCAP     = RHOCI/RHOICE  ! Ice heat capacity (J K-1 Kg-1)
ZEPSILON   = 10._JPRB*EPSILON(ZEPSILON)
ZWHCAP     = 4180._JPRB ! J K-1 Kg-1
ZLN10      = LOG(10.0_JPRB)

! Input fields (in srfsn_driver usually)
  
DO JL=KIDIA,KFDIA
! snow fraction and logical no snow
  ZFRSN(JL)=MAX(PFRTI(JL,5)+PFRTI(JL,7),RFRTINY)

  IF (LDLAND(JL)) THEN
    ! Snowfall is all redirected here ! 
    ZSNOWF(JL)  = PSSFC(JL) + PSSFL(JL)
    ZSNOWF5(JL) = PSSFC5(JL) + PSSFL5(JL)
  ELSE 
    ZSNOWF(JL)  = 0.0_JPRB 
    ZSNOWF5(JL) = 0.0_JPRB 
  ENDIF
! Evap. is only fractional 
  ZEVAPSN(JL) = PFRTI(JL,5)*PEVAPTI(JL,5)  + PFRTI(JL,7)*PEVAPSNW(JL)
  ZEVAPSN5(JL)= PFRTI(JL,5)*PEVAPTI5(JL,5) + PFRTI(JL,7)*PEVAPSNW5(JL)
    
! Rainfall is only fractional 
  ZRAINF(JL) = ZFRSN(JL)*(PTSFC(JL)+PTSFL(JL))
  ZRAINF5(JL)= ZFRSN(JL)*(PTSFC5(JL)+PTSFL5(JL))
  IF(.NOT. YDSOIL%LESNICE)THEN
      ZTBOTTOM5(JL)=PTSURF5(JL)
      ZTBOTTOM(JL) =PTSURF(JL)
  ELSE
    IF (LDLAND(JL))THEN
      ZTBOTTOM5(JL)=PTSURF5(JL)
      ZTBOTTOM(JL) =PTSURF(JL)
    ELSE
      ZTBOTTOM5(JL)=PTICE5(JL)
      ZTBOTTOM(JL) =PTICE(JL)
    ENDIF
  ENDIF

ENDDO


ZZLIQF5=0._JPRB
DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN 
    PTSN(JL,:) = 0._JPRB
    PTSN5(JL,:)= RTT
    PGSN(JL)   = 0.0_JPRB
    PGSN5(JL)  = 0.0_JPRB
    
  ELSE
! heat flux into the snowpack (not scaled by fraction, it simplifies fraction in
! heat transfer eq.
    ZHFLUX(JL)  = PFRTI(JL,5)*(PAHFSTI (JL,5)+YDCST%RLSTT*PEVAPTI (JL,5))&
     & + PFRTI(JL,7)*(PAHFSTI (JL,7)+YDCST%RLVTT*(PEVAPTI (JL,7)-PEVAPSNW (JL))&
     & + YDCST%RLSTT*PEVAPSNW (JL))&
     & + PFRTI(JL,5)*PSSRFLTI (JL,5)&
     & + PFRTI(JL,7)*PSSRFLTI (JL,7)&
     & + (PFRTI(JL,5)+PFRTI(JL,7))*PSLRFL (JL)

    ZHFLUX5(JL) =( PFRTI(JL,5)*(PAHFSTI5(JL,5)+YDCST%RLSTT*PEVAPTI5(JL,5))&
     & + PFRTI(JL,7)*(PAHFSTI5(JL,7)+YDCST%RLVTT*(PEVAPTI5(JL,7)-PEVAPSNW5(JL))&
     & + YDCST%RLSTT*PEVAPSNW5(JL))&
     & + PFRTI(JL,5)*PSSRFLTI5(JL,5)&
     & + PFRTI(JL,7)*PSSRFLTI5(JL,7)&
     & + (PFRTI(JL,5)+PFRTI(JL,7))*PSLRFL5(JL)&
     & )

!*- Soil heat conductivity (simplified) - copied from srfts
    IF (.NOT. YDSOIL%LESNICE) THEN
      ZFF = 0.0_JPRB
      ZWU  = PWSURF(JL)*(1.0_JPRB-ZFF)
      ZWU5 = PWSURF5(JL)*(1.0_JPRB-ZFF)
      ZLWT5 = RLAMBDAWAT**ZWU5
      ZLWT  = ZLWT5*LOG(RLAMBDAWAT)*ZWU

    ZINVWSAT = 1.0_JPRB/RWSATM3D(JL,1)
    ZLIC5 = RLAMBDAICE**(RWSATM3D(JL,1)-ZWU5)
      ZLIC  = -ZLIC5*LOG(RLAMBDAICE)*ZWU
    ZLAMBDASAT  = RLAMSAT1M3D(JL,1)*(ZLIC5*ZLWT+ZLWT5*ZLIC)
    ZLAMBDASAT5 = RLAMSAT1M3D(JL,1)*ZLIC5*ZLWT5
      IF (PWSURF5(JL)*ZINVWSAT > RKERST1) THEN
        ZCOND  = PWSURF (JL)*ZINVWSAT
        ZCOND5 = PWSURF5(JL)*ZINVWSAT
      ELSE
        ZCOND  = 0.0_JPRB
        ZCOND5 = RKERST1
      ENDIF
      ZKERSTEN  = RKERST2*ZCOND/(ZLN10*ZCOND5)
      ZKERSTEN5 = RKERST2*LOG10(ZCOND5)+RKERST3
      ZSOILCOND = (ZLAMBDASAT5-RLAMBDADRYM3D(JL,1))*ZKERSTEN + &
                & ZKERSTEN5*ZLAMBDASAT
      ZSOILCOND5= RLAMBDADRYM3D(JL,1)+ZKERSTEN5*(ZLAMBDASAT5-RLAMBDADRYM3D(JL,1))
      ZSURFCOND = (1._JPRB - PCIL(JL))*ZSOILCOND 
      ZSURFCOND5= (1._JPRB - PCIL(JL))*ZSOILCOND5 + PCIL(JL)*RCONDSICE 
      ZSOILDEPTH1= YDSOIL%RDAW(1) ! 1sth Soil layer depth
    ELSE
      IF (LDLAND(JL))THEN
        ZFF = 0.0_JPRB
        ZWU  = PWSURF(JL)*(1.0_JPRB-ZFF)
        ZWU5 = PWSURF5(JL)*(1.0_JPRB-ZFF)
        ZLWT5 = RLAMBDAWAT**ZWU5
        ZLWT  = ZLWT5*LOG(RLAMBDAWAT)*ZWU

        ZINVWSAT = 1.0_JPRB/RWSATM3D(JL,1)
        ZLIC5 = RLAMBDAICE**(RWSATM3D(JL,1)-ZWU5)
        ZLIC  = -ZLIC5*LOG(RLAMBDAICE)*ZWU
        ZLAMBDASAT  = RLAMSAT1M3D(JL,1)*(ZLIC5*ZLWT+ZLWT5*ZLIC)
        ZLAMBDASAT5 = RLAMSAT1M3D(JL,1)*ZLIC5*ZLWT5
        IF (PWSURF5(JL)*ZINVWSAT > RKERST1) THEN
          ZCOND  = PWSURF (JL)*ZINVWSAT
          ZCOND5 = PWSURF5(JL)*ZINVWSAT
        ELSE
          ZCOND  = 0.0_JPRB
          ZCOND5 = RKERST1
        ENDIF
        ZKERSTEN  = RKERST2*ZCOND/(ZLN10*ZCOND5)
        ZKERSTEN5 = RKERST2*LOG10(ZCOND5)+RKERST3
        ZSOILCOND = (ZLAMBDASAT5-RLAMBDADRYM3D(JL,1))*ZKERSTEN + &
                  & ZKERSTEN5*ZLAMBDASAT
        ZSOILCOND5= RLAMBDADRYM3D(JL,1)+ZKERSTEN5*(ZLAMBDASAT5-RLAMBDADRYM3D(JL,1))
        ZSURFCOND = (1._JPRB - PCIL(JL))*ZSOILCOND 
        ZSURFCOND5= (1._JPRB - PCIL(JL))*ZSOILCOND5 + PCIL(JL)*RCONDSICE 
        ZSOILDEPTH1= YDSOIL%RDAW(1) ! 1sth Soil layer depth
      ELSE
        ZSURFCOND =0._JPRB
        ZSURFCOND5=RCONDSICE
        ZSOILDEPTH1=YDSOIL%RDAI(1)
      ENDIF
     ENDIF
!*-

!! Preparation
    DO JK=1,KLEVSN
      IF (PSSNM1M5(JL,JK) > ZEPSILON ) KLACT=JK
    ENDDO
    DO JK=1,KLEVSN
      ZICE(JK)  = PSSNM1M(JL,JK) - PWSNM1M(JL,JK) ! Ice (Kg m-2)
      ZICE5(JK) = PSSNM1M5(JL,JK) - PWSNM1M5(JL,JK) 

      ! Limit thermally active snow depth to 1 meter of snow (for glaciers in particular)
      ZDSN(JK)  = (PSSNM1M(JL,JK) / PRSNM1M5(JL,JK)) ! read snow depth (m)
      ZDSN5(JK) = (PSSNM1M5(JL,JK)/ PRSNM1M5(JL,JK)) ! read snow depth (m)
      IF (ZDSN5(JK) < RDSNMAX) THEN
        ZDSN(JK)    = (PSSNM1M(JL,JK) / PRSNM1M5(JL,JK))
        ZDSN5(JK)   = PSSNM1M5(JL,JK)/ PRSNM1M5(JL,JK)
      ELSE
        ZDSN(JK)    = 0._JPRB
        ZDSN5(JK)   = RDSNMAX
      ENDIF
      ZTMPHC=(ZICE5(JK)/PRSNM1M5(JL,JK)) 
      IF (ZTMPHC < RDSNMAX) THEN
        ZSNHC(JK)   = (ZIHCAP*PRSNM1M5(JL,JK) * (ZICE(JK)/PRSNM1M5(JL,JK)) + ZWHCAP*PWSNM1M(JL,JK) ) * ZTMST
        ZSNHC5(JK)  = (ZIHCAP*PRSNM1M5(JL,JK) * (ZICE5(JK)/PRSNM1M5(JL,JK))+ ZWHCAP*PWSNM1M5(JL,JK) ) * ZTMST
      ELSE
        ZSNHC(JK)   = ZWHCAP*PWSNM1M(JL,JK) * ZTMST
        ZSNHC5(JK)  = (ZIHCAP*PRSNM1M5(JL,JK)*RDSNMAX + ZWHCAP*PWSNM1M5(JL,JK) ) * ZTMST
      ENDIF

      ! heat conductivity from water vapor transport into the snowpack
      ZTMPCOND=(SNHCONDAV-SNHCONDBV/(PTSNM1M5(JL,JK)-SNHCONDCV)) 
      IF (ZTMPCOND > 0._JPRB) THEN
        ZSNVCOND=SNHCONDPOV*( SNHCONDBV*PAPRS5(JL)*PTSNM1M(JL,JK) - &
                          &(SNHCONDCV-PTSNM1M5(JL,JK))*PAPRS(JL)*&
                          & (SNHCONDAV*SNHCONDCV-SNHCONDAV*PTSNM1M5(JL,JK) + SNHCONDBV) )/&
                          & (PAPRS5(JL)**2._JPRB*(PTSNM1M5(JL,JK)-SNHCONDCV)**2._JPRB)
        ZSNVCOND5=(SNHCONDPOV/PAPRS5(JL))*(SNHCONDAV-SNHCONDBV/(PTSNM1M5(JL,JK)-SNHCONDCV))
      ELSE
        ZSNVCOND5=0._JPRB
        ZSNVCOND =0._JPRB
      ENDIF

      ! snow heat conductivity 
      ZSNCOND(JK)  = ZSNVCOND                           
      ZSNCOND5(JK) = FSNTCOND(PRSNM1M5(JL,JK))+ZSNVCOND5 ! add the thermal cond from water vapor transfer
    ENDDO


   ! Snow conductivity on half levels:
    ZSNCONDH(1)= 4._JPRB*ZFRSN(JL)*(ZDSN5(1)*ZSNCOND(1)-2._JPRB*ZSNCOND5(1)*ZDSN(1))/&
               & ZDSN5(1)**3._JPRB   ! ACTUALY NOT USED ! (W m-2 K-1)
    ZSNCONDH5(1)=ZFRSN(JL)*ZSNCOND5(1) / (0.5_JPRB*ZDSN5(1))*(0.5_JPRB*ZDSN5(1))
    DO JK=2,KLACT
      ZSNCONDH1   = ZFRSN(JL)*2._JPRB*( &
                  & (ZDSN(JK-1)*ZSNCOND5(JK-1)+ZDSN5(JK-1)*ZSNCOND(JK-1))*(ZDSN5(JK-1)+ZDSN5(JK))-&
                  & 2._JPRB*ZDSN5(JK-1)*ZSNCOND5(JK-1)*(ZDSN(JK-1)+ZDSN(JK)) & 
                  & )/(ZDSN5(JK-1)+ZDSN5(JK))**3._JPRB
      ZSNCONDH2   = ZFRSN(JL)*2._JPRB*( &
                  & (ZDSN(JK)*ZSNCOND5(JK)+ZDSN5(JK)*ZSNCOND(JK))*(ZDSN5(JK-1)+ZDSN5(JK))-&
                  & 2._JPRB*ZDSN5(JK)*ZSNCOND5(JK)*(ZDSN(JK-1)+ZDSN(JK)) &
                  & )/(ZDSN5(JK-1)+ZDSN5(JK))**3._JPRB

      ZSNCONDH(JK)= ZSNCONDH1 + ZSNCONDH2
      ZSNCONDH5(JK) = ZFRSN(JL)*2._JPRB*(ZDSN5(JK-1)*ZSNCOND5(JK-1)+ZDSN5(JK)*ZSNCOND5(JK))/&
                & ((ZDSN5(JK-1)+ZDSN5(JK))*(ZDSN5(JK-1)+ZDSN5(JK)))

     ENDDO
    ! Simplified basal heat conduction
     ZSNCONDH(KLACT+1) = ZFRSN(JL)*( ((ZDSN(KLACT)*ZSNCOND5(KLACT)+ZDSN5(KLACT)*ZSNCOND(KLACT)+ZSOILDEPTH1*ZSURFCOND)*(ZDSN5(KLACT)+ZSOILDEPTH1) ) - &
         & (2._JPRB*ZDSN(KLACT)*(ZDSN5(KLACT)*ZSNCOND5(KLACT)+ZSOILDEPTH1*ZSURFCOND5)) )/&
         & (ZDSN5(KLACT)+ZSOILDEPTH1)**3._JPRB

     ZSNCONDH5(KLACT+1)=ZFRSN(JL)*(ZDSN5(KLACT)*ZSNCOND5(KLACT) +  ZSOILDEPTH1*ZSURFCOND5)/&
                  &(ZDSN5(KLACT)+ZSOILDEPTH1)**2._JPRB

    ! Compute coefficients
    IF (KLACT > 1 ) THEN 
      JK=1
      ZTA(JK)=0._JPRB
      ZTA5(JK)=0._JPRB
      ZTB(JK)=ZSNHC(JK)+ZSNCONDH(JK+1)
      ZTB5(JK)=ZSNHC5(JK)+ZSNCONDH5(JK+1)
      ZTC(JK)=-ZSNCONDH(JK+1)
      ZTC5(JK)=-ZSNCONDH5(JK+1)
      ZTR(JK)=(ZHFLUX(JL)-PSNOTRS(JL,JK))+&
             & ZSNHC5(JK)*PTSNM1M(JL,JK) + ZSNHC(JK)*PTSNM1M5(JL,JK)
      ZTR5(JK)=ZHFLUX5(JL)-PSNOTRS5(JL,JK)+ZSNHC5(JK)*PTSNM1M5(JL,JK)
      DO JK=2,KLACT-1
        ZTA(JK) = -ZSNCONDH(JK)
        ZTA5(JK)= -ZSNCONDH5(JK)
        ZTB(JK) = ZSNHC(JK)+ZSNCONDH(JK)+ZSNCONDH(JK+1)
        ZTB5(JK)= ZSNHC5(JK)+ZSNCONDH5(JK)+ZSNCONDH5(JK+1)
        ZTC(JK) = -ZSNCONDH(JK+1)
        ZTC5(JK)= -ZSNCONDH5(JK+1)
        ZTR(JK) = (ZSNHC(JK)*PTSNM1M5(JL,JK) + ZSNHC5(JK)*PTSNM1M(JL,JK)) + &
                & PSNOTRS(JL,JK)
        ZTR5(JK)= ZSNHC5(JK)*PTSNM1M5(JL,JK) + PSNOTRS5(JL,JK)

      ENDDO
! Add radiative flux to second to bottom layer:
      JK=KLACT
      ZTA(JK)  = -ZSNCONDH(JK)
      ZTA5(JK) = -ZSNCONDH5(JK)
      ZTB(JK)  = ZSNHC(JK)+ZSNCONDH(JK)+ZSNCONDH(JK+1)
      ZTB5(JK) = ZSNHC5(JK)+ZSNCONDH5(JK)+ZSNCONDH5(JK+1)
      ZTC(JK)  = 0._JPRB
      ZTC5(JK) = 0._JPRB
      IF(.NOT. YDSOIL%LESNICE)THEN
      ZTR(JK)  = (ZSNHC(JK)*PTSNM1M5(JL,JK)+ZSNHC5(JK)*PTSNM1M(JL,JK))+&
               & (ZSNCONDH(JK+1)*PTSURF5(JL)+ZSNCONDH5(JK+1)*PTSURF(JL))+&
               & PSNOTRS(JL,JK)
      ZTR5(JK) = ZSNHC5(JK)*PTSNM1M5(JL,JK)+ZSNCONDH5(JK+1)*PTSURF5(JL)+PSNOTRS5(JL,JK)
      ELSE
      ZTR(JK)  = (ZSNHC(JK)*PTSNM1M5(JL,JK)+ZSNHC5(JK)*PTSNM1M(JL,JK))+&
               & (ZSNCONDH(JK+1)*ZTBOTTOM5(JL)+ZSNCONDH5(JK+1)*ZTBOTTOM(JL))+&
               & PSNOTRS(JL,JK)
      ZTR5(JK) = ZSNHC5(JK)*PTSNM1M5(JL,JK)+ZSNCONDH5(JK+1)*ZTBOTTOM5(JL)+PSNOTRS5(JL,JK)
      ENDIF
    ENDIF
    
! SOLVE
    ZTSTAR(:)  = PTSNM1M(JL,:)
    ZTSTAR5(:) = PTSNM1M5(JL,:)
    ZGSN = 0._JPRB
    ZGSN5= 0._JPRB
    IF (KLACT == 1 ) THEN
      ! SINGLE LAYER
      ZTB(1) = ZSNHC(1)+ZSNCONDH(2)
      ZTB5(1)= ZSNHC5(1)+ZSNCONDH5(2)
      IF(.NOT. YDSOIL%LESNICE)THEN
      ZTR(1) = (ZHFLUX(JL)-PSNOTRS(JL,1))+&
             & ZSNHC(1)*PTSNM1M5(JL,1)+ZSNHC5(1)*PTSNM1M(JL,1)+&
             & ZSNCONDH(2)*PTSURF5(JL)+ZSNCONDH5(2)*PTSURF(JL)
      ZTR5(1)= (ZHFLUX5(JL)-PSNOTRS5(JL,1))+ZSNHC5(1)*PTSNM1M5(JL,1)+ZSNCONDH5(2)*PTSURF5(JL)
      ELSE
      ZTR(1) = (ZHFLUX(JL)-PSNOTRS(JL,1))+&
             & ZSNHC(1)*PTSNM1M5(JL,1)+ZSNHC5(1)*PTSNM1M(JL,1)+&
             & ZSNCONDH(2)*ZTBOTTOM5(JL)+ZSNCONDH5(2)*ZTBOTTOM(JL)
      ZTR5(1)= (ZHFLUX5(JL)-PSNOTRS5(JL,1))+ZSNHC5(1)*PTSNM1M5(JL,1)+ZSNCONDH5(2)*ZTBOTTOM5(JL)
      ENDIF


      ZTSTAR(1) = (ZTR(1)*ZTB5(1)-ZTB(1)*ZTR5(1))/ZTB5(1)**2._JPRB
      ZTSTAR5(1)=ZTR5(1)/ZTB5(1)
      IF(.NOT. YDSOIL%LESNICE)THEN
      ZGSN = ZSNCONDH(2)*(ZTSTAR5(1) -PTSURF5(JL))+ZSNCONDH5(2)*(ZTSTAR(1)-PTSURF(JL))
      ZGSN5= ZSNCONDH5(2)*(ZTSTAR5(1)-PTSURF5(JL))
      ELSE
      ZGSN = ZSNCONDH(2)*(ZTSTAR5(1) -ZTBOTTOM5(JL))+ZSNCONDH5(2)*(ZTSTAR(1)-ZTBOTTOM(JL))
      ZGSN5= ZSNCONDH5(2)*(ZTSTAR5(1)-ZTBOTTOM5(JL))
      ENDIF
    ELSE
    ! MULTI LAYER
      CALL TRISOLVERSTL(KLACT,ZTA5,ZTB5,ZTC5,ZTR5,ZTSTAR5,&
                       &ZTA,ZTB,ZTC,ZTR,ZTSTAR)
      IF(.NOT. YDSOIL%LESNICE)THEN
      ZGSN = ZSNCONDH(KLACT+1)*ZTSTAR5(KLACT) + ZSNCONDH5(KLACT+1)*ZTSTAR(KLACT) - &
           & (ZSNCONDH(KLACT+1)*PTSURF5(JL)+ZSNCONDH5(KLACT+1)*PTSURF(JL))
      ZGSN5=ZSNCONDH5(KLACT+1)*(ZTSTAR5(KLACT)-PTSURF5(JL))
      ELSE
      ZGSN = ZSNCONDH(KLACT+1)*ZTSTAR5(KLACT) + ZSNCONDH5(KLACT+1)*ZTSTAR(KLACT) - &
           & (ZSNCONDH(KLACT+1)*ZTBOTTOM5(JL)+ZSNCONDH5(KLACT+1)*ZTBOTTOM(JL))
      ZGSN5=ZSNCONDH5(KLACT+1)*(ZTSTAR5(KLACT)-ZTBOTTOM5(JL))
      ENDIF
    ENDIF 

    ! * SPECIAL CASE WHEN THERE IS MELTING ON THE 1ST LAYER
    IF ( ZTSTAR5(1) > RTT ) THEN 
      IF (KLACT == 1 ) THEN
        ! SINGLE LAYER
        IF(.NOT. YDSOIL%LESNICE)THEN
        ZGSN  = ZSNCONDH(2)*(RTT-PTSURF5(JL)) - ZSNCONDH5(2)*PTSURF(JL)
        ZGSN5 = ZSNCONDH5(2)*(RTT-PTSURF5(JL))
        ELSE
        ZGSN  = ZSNCONDH(2)*(RTT-ZTBOTTOM5(JL)) - ZSNCONDH5(2)*ZTBOTTOM(JL)
        ZGSN5 = ZSNCONDH5(2)*(RTT-ZTBOTTOM5(JL))
        ENDIF
        ZTSTAR(1) = -ZSNHC(1)*ZHFLUX5(JL)/ZSNHC5(1)**2._JPRB + ZHFLUX(JL)/ZSNHC5(1) - &
                  & (-ZSNHC(1)*PSNOTRS5(JL,1)/ZSNHC5(1)**2._JPRB +PSNOTRS(JL,1)/ZSNHC5(1)) + &
                  & PTSNM1M(JL,1) - &
                  & (-ZSNHC(1)*ZGSN5/ZSNHC5(1)**2._JPRB + ZGSN/ZSNHC5(1))
        ZTSTAR5(1)= (ZHFLUX5(JL)-PSNOTRS5(JL,1)+ZSNHC5(1)*PTSNM1M5(JL,1)-ZGSN5)/ZSNHC5(1)
      ELSE
        ! MULTI LAYER
        ! SOLVE SYSTEM SETTING T1 == RTT 
        ZTB (1)  = 0._JPRB
        ZTB5(1) = 1._JPRB
        ZTC (1)  = 0._JPRB 
        ZTC5(1) = 0._JPRB 
        ZTR (1)  = 0._JPRB
        ZTR5(1) = RTT 
        ZTSTAR(:)  = PTSNM1M(JL,:)  
        ZTSTAR5(:) = PTSNM1M5(JL,:)  
        CALL TRISOLVERSTL(KLACT,ZTA5,ZTB5,ZTC5,ZTR5,ZTSTAR5,&
                         &ZTA,ZTB,ZTC,ZTR,ZTSTAR)
        IF(.NOT. YDSOIL%LESNICE)THEN
        ZGSN = ZSNCONDH(KLACT+1)*(ZTSTAR5(KLACT)-PTSURF5(JL))+ZSNCONDH5(KLACT+1)*(ZTSTAR(KLACT)-PTSURF(JL))
        ZGSN5= ZSNCONDH5(KLACT+1)*(ZTSTAR5(KLACT)-PTSURF5(JL))
        ELSE
        ZGSN = ZSNCONDH(KLACT+1)*(ZTSTAR5(KLACT)-ZTBOTTOM5(JL))+ZSNCONDH5(KLACT+1)*(ZTSTAR(KLACT)-ZTBOTTOM(JL))
        ZGSN5= ZSNCONDH5(KLACT+1)*(ZTSTAR5(KLACT)-ZTBOTTOM5(JL))
        ENDIF
        ! EXPLICIT 1ST LAYER TEMP CALCULATION
        ZTSTAR(1) = -ZSNHC(1)*ZHFLUX5(JL)/ZSNHC5(1)**2._JPRB + ZHFLUX(JL)/ZSNHC5(1) - &
                  & (-ZSNHC(1)*PSNOTRS5(JL,1)/ZSNHC5(1)**2._JPRB +PSNOTRS(JL,1)/ZSNHC5(1)) + &
                  & PTSNM1M(JL,1) - &
                  & (ZSNCONDH(2)*ZTSTAR5(1) + ZSNCONDH5(2)*ZTSTAR(1) - ZSNCONDH(2)*ZTSTAR5(2)- ZSNCONDH5(2)*ZTSTAR(2) )
        ZTSTAR5(1) = (ZHFLUX5(JL)-PSNOTRS5(JL,1)+ZSNHC5(1)*PTSNM1M5(JL,1)-ZSNCONDH5(2)*(ZTSTAR5(1)-ZTSTAR5(2)))/ZSNHC5(1)
      ENDIF
    ENDIF
    
!! UPDATE SOLID/LIQUID CONTENTS OF FIRST LAYER BEFORE FURTHER CALCULATIONS
    ZISTAR(:)  = ZICE(:)
    ZISTAR5(:) = ZICE5(:)
    ZWSTAR(:)  = PWSNM1M(JL,:)
    ZWSTAR5(:) = PWSNM1M5(JL,:)
    ZLIQF(:)   = 0._JPRB
    ZLIQF5(:)  = 0._JPRB
      !! INCLUDE SNOWFALL AND RAINFALL INTERCEPTION 
    ZISTAR(1)  = ZICE(1)+PTMST*(ZSNOWF(JL)+ZEVAPSN(JL))
    ZISTAR5(1) = ZICE5(1)+PTMST*(ZSNOWF5(JL)+ZEVAPSN5(JL))
    ZLIQF(0)   = PTMST*ZRAINF(JL)
    ZLIQF5(0)  = PTMST*ZRAINF5(JL)
    
    !! PHASE CHANGES AND LIQUID WATER BALANCE STARTING IN LAYER 1
       ZMELTSN(1:KLEVSN)  = 0._JPRB
       ZMELTSN5(1:KLEVSN) = 0._JPRB
       ZFREZSN(1:KLEVSN)  = 0._JPRB
       ZFREZSN5(1:KLEVSN) = 0._JPRB
       ZGSNRES(0:KLEVSN)  = 0._JPRB
       ZGSNRES5(0:KLEVSN) = 0._JPRB
       ZQ(1:KLEVSN)       = 0._JPRB
       ZQ5(1:KLEVSN)      = 0._JPRB
       ZTMP05(1:KLEVSN)   = 0._JPRB
       ZTMP0(1:KLEVSN)   = 0._JPRB
    DO JK=1,1
      ! UPDATE LIQUID WATER FROM FLUX ABOVER 
      ZWSTAR(JK)  = PWSNM1M(JL,JK)  + ZLIQF(JK-1)
      ZWSTAR5(JK) = PWSNM1M5(JL,JK) + ZLIQF5(JK-1)
      
      ! PHASE CHANGES 
      ZTMP0(JK)   = ZSNHC(JK)*ZTSTAR5(JK)+ZSNHC5(JK)*ZTSTAR(JK)-ZSNHC(JK)*RTT
      ZTMP05(JK)  = ZSNHC5(JK)*(ZTSTAR5(JK)-RTT)
      ZCHANGE5= RLMLT*ZTMST*ZISTAR5(JK)
      IF (ZTMP05(JK) < ZCHANGE5 ) THEN
         ZMELTSN(JK)=ZTMP0(JK)
         ZMELTSN5(JK)=ZTMP05(JK)
      ELSE
         ZMELTSN(JK) =RLMLT*ZTMST*ZISTAR(JK)
         ZMELTSN5(JK)=RLMLT*ZTMST*ZISTAR5(JK)
      ENDIF
      IF (ZMELTSN5(JK)>0._JPRB) THEN
         ZMELTSN(JK)=ZMELTSN(JK)
         ZMELTSN5(JK)=ZMELTSN5(JK)
      ELSE
         ZMELTSN(JK)=0._JPRB
         ZMELTSN5(JK)=0._JPRB
      ENDIF

      ZCHANGE5 = -RLMLT*ZTMST*ZWSTAR5(JK)
      IF (ZTMP05(JK) > ZCHANGE5) THEN
        ZFREZSN(JK)  = ZTMP0(JK)
        ZFREZSN5(JK) = ZTMP05(JK)
      ELSE
        ZFREZSN(JK)  = -RLMLT*ZTMST*ZWSTAR(JK)
        ZFREZSN5(JK) = -RLMLT*ZTMST*ZWSTAR5(JK)
      ENDIF
      IF (ZFREZSN5(JK) > 0._JPRB) THEN
        ZFREZSN(JK)  = 0._JPRB
        ZFREZSN5(JK) = 0._JPRB
      ELSE 
        ZFREZSN(JK)  = ZFREZSN(JK)
        ZFREZSN5(JK) = ZFREZSN5(JK)
      ENDIF

      ZQ(JK)    = ZMELTSN(JK) + ZFREZSN(JK)
      ZQ5(JK)   = ZMELTSN5(JK) + ZFREZSN5(JK)
      
      ! FINAL TEMP. UPDATE
      ZTSTAR(JK)  = ZTSTAR(JK)  - &
                  & ((ZQ(JK)+ZGSNRES(JK-1))*ZSNHC5(JK)-ZSNHC(JK)*(ZQ5(JK)+ZGSNRES5(JK-1)))/ZSNHC5(JK)**2._JPRB
      ZTSTAR5(JK) = ZTSTAR5(JK) - (ZQ5(JK)+ZGSNRES5(JK-1))/ZSNHC5(JK)

      IF ( ZTSTAR5(JK) > RTT ) THEN
        ZGSNRES(JK) = ZSNHC(JK)*RTT - &
                & (ZSNHC(JK)*ZTSTAR5(JK) + ZSNHC5(JK)*ZTSTAR(JK))
        ZGSNRES5=ZSNHC5(JK)*(RTT-ZTSTAR5(JK))
        ZTSTAR(JK) =0._JPRB
        ZTSTAR5(JK)=RTT
      ELSE
        ZGSNRES(JK) =0._JPRB
        ZGSNRES5(JK)=0._JPRB
      ENDIF
      
      ! UPDATE SOLID / LIQUID MASS
      ZWSTAR(JK)  = ZWSTAR(JK)  + ZQ(JK)/RLMLT*PTMST 
      ZWSTAR5(JK) = ZWSTAR5(JK) + ZQ5(JK)/RLMLT*PTMST 
      ZISTAR(JK)  = ZISTAR(JK)  - ZQ(JK)/RLMLT*PTMST
      ZISTAR5(JK) = ZISTAR5(JK) - ZQ5(JK)/RLMLT*PTMST
      
      ! LIQUD WATER BUDGET
      ! SNOW LIQUID WATER CAPACITY
      ZCONSLWC=MAX(0._JPRB,(RLWCSWEA-PRSNM1M5(JL,JK))/RLWCSWEA)
      ZWCAP   = ZISTAR(JK)*(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC)
      ZWCAP5  = ZISTAR5(JK)*(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC)
  
      ZLIQF5(JK) = (ZWSTAR5(JK)-ZWCAP5)
      IF (ZLIQF5(JK)>0._JPRB) THEN
        ZLIQF(JK)  = (ZWSTAR(JK)-ZWCAP)
        ZLIQF5(JK) = (ZWSTAR5(JK)-ZWCAP5)
      ELSE
        ZLIQF(JK)  = 0._JPRB
        ZLIQF5(JK) = 0._JPRB
      ENDIF

      ENDDO

      ! Fill unused levels
      DO JK=KLACT+1,KLEVSN
       ! Original used in tl tests
       ! ZTSTAR(JK) = ZTSTAR(KLACT)
        ZTSTAR(JK) = 0._JPRB
        ZTSTAR5(JK)= ZTSTAR5(KLACT)
      ENDDO  
      
    !! FINAL VALUES
    DO JK=1,KLEVSN
      IF (ZTSTAR5(JK)<RTT) THEN
        PTSN(JL,JK) = ZTSTAR(JK)
        PTSN5(JL,JK)= ZTSTAR5(JK)
      ELSE
        PTSN(JL,JK) = 0._JPRB
        PTSN5(JL,JK)= RTT
      ENDIF
    ENDDO
    PGSN(JL) = ZGSN - ZGSNRES(KLACT) + PSNOTRS(JL,KLACT+1)
    PGSN5(JL)= ZGSN5- ZGSNRES5(KLACT)+ PSNOTRS5(JL,KLACT+1)

    DO JK=1,KLEVSN
      IF ((PTSN5(JL,JK)<100._JPRB)) THEN
        write(NULERR,*) 'Very cold snow temperature, webalstl'
        write(NULERR,*) 'Tsn-1',PTSNM1M5(JL,:)
        write(NULERR,*) 'Tsn',PTSN5(JL,:)
        write(NULERR,*) 'SWE-1',PSSNM1M5(JL,:)
        PTSN(JL,JK)=0._JPRB
        PTSN5(JL,JK)=100.0_JPRB
        !*CALL ABORT_SURF('Very snow cold temperature')
      ENDIF
    ENDDO 

  ENDIF 
ENDDO
  
END ASSOCIATE
!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALSTL_MOD:SRFSN_WEBALSTL',1,ZHOOK_HANDLE)

END SUBROUTINE SRFSN_WEBALSTL

SUBROUTINE TRISOLVERSTL(KSNLAC,ZA5,ZB5,ZC5,ZF5,ZTSTAR5,&
                       &ZA,ZB,ZC,ZF,ZTSTAR)
USE PARKIND1, ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
IMPLICIT NONE

!! SINGLE POINT TRIDIAGONAL SOLVER ! 

!* DECLARATION OF ARGUMENTS
INTEGER(KIND=JPIM)  ,INTENT(IN)  :: KSNLAC ! NUMBER OF ACTIVE SNOW LAYERS
REAL(KIND=JPRB),INTENT(IN),DIMENSION(:)  :: ZA5,ZB5,ZC5,ZF5
REAL(KIND=JPRB),INTENT(IN),DIMENSION(:)  :: ZA, ZB, ZC, ZF

REAL(KIND=JPRB),INTENT(INOUT),DIMENSION(:) :: ZTSTAR5
REAL(KIND=JPRB),INTENT(INOUT),DIMENSION(:) :: ZTSTAR

!* LOCAL VARIABLES  
INTEGER(KIND=JPIM)  :: JK
REAL(KIND=JPRB),DIMENSION(SIZE(ZB)) :: BET5,BET
REAL(KIND=JPRB),DIMENSION(SIZE(ZB)) :: GAM5,GAM

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALSTL_MOD:TRISOLVERSTL',0,ZHOOK_HANDLE)

! top layer elimination
BET(1)=ZB(1)
BET5(1)=ZB5(1)

ZTSTAR(1) = (ZF(1)*BET5(1) - ZF5(1)*BET(1))/BET5(1)**2._JPRB
ZTSTAR5(1)= ZF5(1)/BET5(1)

! First downward scan
DO JK=2,KSNLAC
  GAM5(JK)= ZC5(JK-1)/BET5(JK-1)
  GAM(JK) = (ZC(JK-1)*BET5(JK-1) - ZC5(JK-1)*BET(JK-1))/BET5(JK-1)**2._JPRB

  BET5(JK)= ZB5(JK)-ZA5(JK)*GAM5(JK)
  BET(JK) = ZB(JK)-(ZA(JK)*GAM5(JK) + ZA5(JK)*GAM(JK))

  ZTSTAR(JK) = (ZF(JK)*BET5(JK) - ZF5(JK)*BET(JK))/BET5(JK)**2._JPRB - &
             & (BET5(JK)*(ZA(JK)*ZTSTAR5(JK-1)+ZA5(JK)*ZTSTAR(JK-1))-BET(JK)*(ZA5(JK)*ZTSTAR5(JK-1)))/BET5(JK)**2._JPRB
  ZTSTAR5(JK)= (ZF5(JK)-ZA5(JK)*ZTSTAR5(JK-1))/BET5(JK)
ENDDO
! Back substitution
DO JK=KSNLAC-1,1,-1
  ZTSTAR(JK)  = ZTSTAR(JK)-&
              & (GAM(JK+1)*ZTSTAR5(JK+1)+GAM5(JK+1)*ZTSTAR(JK+1))
  ZTSTAR5(JK) = ZTSTAR5(JK)-GAM5(JK+1)*ZTSTAR5(JK+1)
ENDDO

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALSTL_MOD:TRISOLVERTL',1,ZHOOK_HANDLE)
END SUBROUTINE TRISOLVERSTL


END MODULE SRFSN_WEBALSTL_MOD
